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PACS 81.10.Aj - Theory and models of crystal growth; physics of crystal growth, crystal mor- 
phology and orientation 
PACS 64 . 70 . Dv - Solid-liquid transitions 
PACS 02 . 70 . Ns - Molecular dynamics and particle methods 

Abstract. - The melting and crystallization of AlsoNiso are studied by means of molecular 
dynamics computer simulations, using a potential of the embedded atom type to model the inter- 
actions between the particles. Systems in a slab geometry are simulated where the B2 phase of 
AINi in the middle of an elongated simulation box is separated by two planar interfaces from the 
liquid phase, thereby considering the (100) crystal orientation. By determining the temperature 
dependence of the interface velocity, an accurate estimate of the melting temperature is provided. 
The value k — 0.0025 m/s/K for the kinetic growth coefficient is found. This value is about two 
orders of magnitude smaller than that found in recent simulation studies of one-component metals. 
The classical Wilson-Frenkel model is not able to describe the crystal growth kinetics on a quan- 
titative level. We argue that this is due to the neglect of diffusion processes in the liquid-crystal 
interface. 
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Introduction. — The classical model for crystalliza- 
tion from the melt is the one proposed by Wilson [1] and 
Frcnkcl [2]. It considers crystal growth as an activated 
process, controled by the mass transport in the liquid. 
However, various studies using molecular dynamics (MD) 
computer simulation have shown that the Wilson-Frenkel 
scenario is not applicable to a large class of materials. 
Especially in pure metals, growth kinetics is much faster 
than expected for an activated diffusion-limited mecha- 
nism [3-14]. In this case, rearrangements in the liquid 
structure are not required to provide the formation of crys- 
talline layers. This may explain why one-component met- 
als are not glassforming systems in general. On the other 
hand, binary metallic alloys are known as glassforming 
systems, provided that heterogeneous nucleation can be 
avoided. Indeed, these systems exhibit in general a much 
slower growth kinetics than pure metals. Recent stud- 
ies have demonstrated that the MD simulation technique 
is well-suited to elucidate the crystallization kinetics in 
binary alloys [15-21]. But these studies show also that 
the growth kinetics in binary mixtures is more compli- 
cated than in the one-component counterparts. One of 



the open questions is to what extent the Wilson-Frenkel 
picture is valid for two-component metals. This question 
is addressed in the following. 

In this work, the crystal growth kinetics of the binary al- 
loy AI50N150 is investigated by MD simulation. The exper- 
imental melting temperature for this system is at 1920 K 
where it exhibits a first order phase transition from a liq- 
uid to an intermetallic B2 phase. Very recently, the crys- 
tal growth velocity for this transition has been measured 
by Reutzel et al. [22] using an electromagnetic levitation 
technique under reduced gravity conditions in combination 
with a high-speed camera. At an undercooling of about 
60 K, growth velocities of the order of 0.1 m/s were found. 
This value is about two orders of magnitude smaller than 
that found for pure metals at comparable undercoolings, 
indicating that AI50M50 may be the prototype of a system 
with a diffusion-limited growth mechanism. 

The MD simulation allows for an accurate determina- 
tion of the melting temperature, kinetic growth coeffi- 
cients, and transport coefficients such as self- and interdif- 
fusion constants. These information are required to check 
the validity of the Wilson-Frenkel model of crystal growth. 
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Fig. 1: Snapshot of a simulated configuration with two crystal- 
melt interfaces of the system AlsoNiso at the temperature T — 
1500 K. Al and Ni atoms are shown as grey and brown spheres, 
respectively. 



As we shall sec below the kinetic growth coefficient, as es- 
timated by our simulation for the (100) orientation of the 
crystal, is indeed much smaller than that found for sim- 
ple metals. Thereby, the growth velocities are in good 
agreement with those measured in the aforementioned ex- 
periment by Reutzel et al. [22]. However, we demonstrate 
that the Wilson-Frcnkcl model is not able to describe the 
crystal growth kinetics in AlsoNiso on a quantitative level, 
at least for plausible choices of the various free parameters 
appearing in the theory. This indicates the need of micro- 
scopic theories on the various aspects of crystal growth 
kinetics. By computing diffusion profiles for an inhomo- 
gencous crystal-liquid system at coexistence, we show ex- 
plicitcly that such theories have to take into account dif- 
fusion processes in the crystal-liquid interface region. 

Details of the simulation. — To investigate the crys- 
tallization of AI50M50 from the melt, we have done ex- 
tensive molecular dynamics computer simulations. The 
interactions between the atoms were modelled by a po- 
tential of the embedded atom type, proposed by Mishin 
et al. [23]. Recent studies have shown that this potential 
gives a realistic description of the diffusion dynamics in 
Al-Ni melts [24,25]. The simulations were done at con- 
stant pressure (p C xt = 0). For this, an algorithm proposed 
by Andersen was used, setting the mass of the piston to 
0.0027 u [26]. Temperature was kept constant by coupling 
the system at every 100 steps to a stochastic heat bath. 
The equations of motion were integrated with the velocity 
form of the Verlet algorithm [27] with a time step of 1 fs. 

At each temperature in the range 1600 K> T > 1200 K, 
12 independent samples with solid- liquid interfaces were 
prepared. To this end, the B2 phase of A1.5oNi 5 o was 
equilibrated at the target temperature for 1 ns. The sim- 
ulations were done for a system of N = 3072 particles 
(A r Ai = Nm = 1536) in an elongated simulation box of size 
L x L x L z (with L z = 3 x L), considering the (100) di- 
rection of the crystal. Periodic boundary conditions were 
employed in all three spatial directions. Having relaxed 
the crystal sample, one third of the particles in the mid- 
dle of the box were fixed and the rest of the system was 
melted during 500 ps at T = 3000 K. Then, the whole sys- 
tem was annealed at the target temperature for another 




Fig. 2: Number density profiles during crystal growth at T = 
1460 K for Al (solid lines) and Ni (dashed lines). The profiles 
corresponding to t = 560 ps and t = 994 ps are shifted with 
respect to the t = Ops profiles by 0.3 A and 0.6 A, as indicated. 



500 ps, before we started the production runs over Ins in 
the NpT ensemble. A snapshot of the system with two 
interfaces at T = 1500 K is shown in Fig. 1. We did also 
20 independent microcanonical runs of the crystal-liquid 
system at the coexistence temperature T = 1520 K (see 
below), starting from fully equilibrated samples. These 
runs, each of them over 1 ns, was used to study the diffu- 
sion dynamics in the crystal-liquid interface region. 

In addition, simulations of liquid samples were per- 
formed at the temperatures T = 1200 K, 1300 K, 1400 K, 
1500 K, 1600 K, 1800 K, and 2000 K, in order to determine 
self-diffusion coefficients as well as the intcrdiffusion coef- 
ficient (see below). In this case, systems of 2000 particles 
were placed in a cubic simulation box. At each tempera- 
ture, equilibration runs over 1 ns were done in the NpT en- 
semble, followed by microcanonical production runs over 
23 ns. 

Results. — As described in the previous section, sam- 
ples in an elongated simulation box were prepared as start- 
ing configurations where the crystal in the middle is sur- 
rounded by the liquid phase on both sides, separated by 
two interfaces (see Fig. 1). The behavior of these sam- 
ples depends strongly on the temperature at which they 
are simulated. While below the melting temperature T m , 
the crystal will grow (as shown in Fig. 2), it will melt 
above T m . From the simulation, the velocity vi with which 
the liquid-crystal interface moves can be determined. At 
T = T m , the interface velocity v\ vanishes. Thus, by the 
extrapolation v\ — > the melting temperature T m can be 
estimated. In the following, we show that this procedure 
yields a rather accurate estimate of T m . Then, we demon- 
strate that the crystal growth mechanism in AI50M50 can 
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Fig. 3: Volume per particle, V p as a function of time for dif- 
ferent temperatures, as indicated. The bold dashed lines are 
examples of linear fits from which the volume velocity V is 
determined. 



be elucidated by investigating the diffusion dynamics in 
the liquid phase and in the crystal-liquid interface region. 

Figure 2 displays the partial number density profiles 
p(z) of Al and Ni at T = 1460 K along the z direction, 
i.e. perpendicular to the solid-liquid interfaces. The lower 
profiles in Fig. 2 correspond to the starting configuration, 
while the second and the third ones correspond to t = 
560 ps and 994 ps. Note that in Fig. 2 the z coordinate is 
scaled by the factor 2/L z , placing z = in the middle of 
the simulation box. Whereas the crystal structure leads to 
pronounced peaks in p(z), a constant density is observed 
for the liquid regions along the z direction, as expected. 
We can also infer from Fig. 2 that the intermetallic B2 
phase [here in (100) orientation] exhibits a pronounced 
chemical ordering, characterized by the alternate sequence 
of Al and Ni layers. This indicates that, different from 
one-component metals, the crystal growth kinetics relies 
on local rearrangements in the liquid structure. Thus, one 
may expect that diffusive transport is required to bring the 
atoms of each species to a suitable site in the B2 crystal. 
As one can further see in Fig. 2, the crystal is growing at 
T = 1460 K. Thus, this temperature is below the melting 
temperature of our AI50M50 model. 

Since the density of the crystalline B2 phase is higher 
than that of the liquid phase, the total volume of the 
system decreases at temperatures T < T m whereas it in- 
creases above T m . Figure 3 shows the time dependence 
of the volume per particle, V p , for different temperatures 
between 1400 K and 1600 K. From this plot, one can in- 
fer that the melting temperature is between 1500 K and 
1530 K. Also shown in Fig. 3 are examples of linear fits 
of the form f(t) = A — V p t. Such linear growth laws are 
expected for steady state growth [7]. We use these fits 
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Fig. 4: Interface velocity as a function of temperature (filled 
circles, the dashed line is a guide to the eye). The bold line is a 
linear fit, vi = fc(T m — T), yielding the kinetic growth coefficient 
k = 0.0025 m/s/K and the melting temperature T m = 1520 K. 
The inset shows the interface velocity as a function of under- 
cooling T m — T and the interface velocity divided by the aver- 
aged self-diffusion coefficient (dashed line). 



to determine the change of the volume V per unit time. 
The deviations from the linear behavior at short times re- 
veal that the growth (or melting) of the crystal is not yet 
in a steady state regime [7]. At high temperatures, we 
see a complete melting of the crystal and thus the volume 
V p reaches a constant at long times corresponding to the 
specific volume of the liquid phase. Prior to this, the melt- 
ing of the crystal is faster than in the linear steady-state 
regime. In this intermediate regime the crystal has shrunk 
to such small dimensions that we see essentially the inter- 
action between the two interfaces in the simulation box 
and thus strong deviations from steady state growth are 
observed. 

From the volume change Vp, the velocity v\, with which 
the liquid-crystal interfaces move, can be estimated as fol- 
lows: 



''1 



2Ni(V c - Vi) 



(1) 



Here, the product N\(V C — V\) quantifies the increase of 
the volume caused by the addition of a crystalline layer 
(with N\ the average number of particles in a layer, and 
V c and Vi the specific volumes of the crystal and the liquid 
phase, respectively). The length d is the spacing between 
crystalline layers. 

Figure 4 displays the interface velocity v\ as a function 
of temperature. We see that v\ vanishes around 1520 K 
and thus this temperature is the estimate for the melt- 
ing temperature, T m , of our simulation model. Note that 
the experimental value for T m is around 1920 K and so 
our simulation underestimates the experimental value by 
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about 20%. Around T m , the simulation data for v\ can be 
fitted by the linear law vi = fc(T m — T) where the fit pa- 
rameter k is the so-called kinetic coefficient. The fit, that 
is shown in Fig. 4, yields the value k = 0.0025 m/s/K. This 
value is about two orders of magnitude smaller than the 
typical values for kinetic coefficients that have been found 
in simulations of one-component metals [9-11,28]. 

The inset in Fig. 4 shows the interface velocity as a 
function of undercooling AT = T m — T. We see that v\ 
increases linearly up to an undercooling of about 30 K. 
At AT Kt 180 K, the interface velocity reaches a maxi- 
mum value of about 0.15 m/s. Note that at small under- 
coolings our simulation data are in good agreement with 
recent experimental data on AI50N150, measured under re- 
duced gravity conditions during a parabolic flight cam- 
paign [22]. Also shown in the inset of Fig. 4 is the quan- 
tity 2i;i/(£>Ni + -Dai), with Z?Ni and Dai the self-diffusion 
constants of Ni and Al, respectively. The self-diffusion 
constants will be discussed in detail below. Here, we note 
that the maximum in v\ disappears when one divides this 
quantity by the averaged self-diffusion coefficient. Thus, 
the occurrence of a maximum in t>i(AT) is due to the 
slowing down of diffusion processes with decreasing tem- 
perature. 

On a qualitative level, the behavior of vi(AT) can be 
understood in the framework of the Wilson-Frenkcl model. 
The model relates the interface velocity to the difference 
between the rate at which the atoms join the crystal and 
the rate at which they leave the crystal. As a result the 
following formula for V\ is obtained [28]: 



vj = A 



kin 



1 — cxp 



Aff 
k R T 



(2) 



with Akin a kinetic prefactor, fee the Boltzmann con- 
stant and Ag the free energy difference between the liquid 
and crystal phase. Close to coexistence, the free energy 
difference Ag is proportional to AT, and the exponen- 
tial function in Eq. (2) can be approximated, such that 
1 — cxp (— 77HO ~ with I the latent heat of the 



liquid-to-solid transition. Furthermore, the kinetic pref- 
actor Akin can be expressed in terms of the diffusion co- 
efficient D of the liquid. Eventually, at small AT the 
expression for vi can be written as [28] 



/cwfAT with 



kwF = -tjD 3) 



where / represents the fraction of collisions with the crys- 
tal that contribute to the growth of the crystal. The pa- 
rameter A corresponds to an elementary diffusive jump 
distance of particles in the liquid [28]. Note that it is 
assumed in the derivation of Eq. (3) that the diffusion 
constant can be expressed by an Arrhenius law, 



D = D exp 



Q 
k B T 



with 
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with Q an activation energy associated with the diffusion 



Fig. 5: Arrhenius plot of self-diffusion constants D a (a = 
Al, Ni), Onsager coefficient L and interdiffusion constant Dab 
for AlsoNiso. The solid lines are fits with the Arrhenius law 
(4), see text. The arrow indicates the location of the melting 
temperature T m ~ 1520 K of the simulation model. 



of the atoms in the liquid and v a frequency of the order 
of the Debye frequency. 

In order to check whether the Wilson-Frenkel formula 
for the kinetic coefficient fcwF m Eq. (3) yields quanti- 
tative agreement with the value k = 0.0025 m/s/K for 
AlsoNiso (see above), we have computed the temperature 
dependence of self- and interdiffusion coefficients around 
T,„ = 1520 K. Whereas the self-diffusion constant D a is 
the transport coefficient for tagged particle diffusion of 
atoms of type a (here a = Al, Ni) , the interdiffusion 
coefficient Dab describes diffusive transport due to con- 
centration fluctuations among the different components. 
The self-diffusion constants D a have been computed from 
the long-time limit of the corresponding mean-squared 
displacements. The interdiffusion coefficient is given by 
Dab = $L where <£> is the so-called thermodynamic fac- 
tor and L is the Onsager coefficient. The thermodynamic 
factor expresses the thermodynamic forces to homogenize 
the mixture with respect to concentration fluctuations. 
We have calculated this quantity from the q — » limit 
of the inverse concentration-concentration structure fac- 
tor (see Ref. [25]). The Onsager coefficient L contains all 
the kinetic contributions to T>ab and can be determined 
from a generalized mean-squared displacement describing 
the centre-of-mass motion of one species. For details of 
the calculation of L and <&, we refer the reader to a recent 
publication [25]. 

An Arrhenius plot of the different diffusion coefficients 
is shown in Fig. 5. As in a recent simulation study of 
Al 80 Ni 2 o [25], the interdiffusion coefficient is about a factor 
4 to 6 higher than the self-diffusion constants. This is due 
to the thermodynamic factor (note that the Onsager coef- 
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Fig. 6: Number density profiles and diffusion profiles for Ni 
and Al at coexistence, as indicated. 



ficient lies below the self-diffusion constants). The origin 
of this behavior is a large resistance to macroscopic con- 
centration fluctuations in dense liquids (a similar property 
of dense liquids is their very low compressibility). With 
respect to crystal growth kinetics in a binary alloy such as 
AI50M50, it is not clear whether one has to consider self- or 
intcrdiffusive transport as the limiting growth mechanism. 

Also shown in Fig. 5 are fits with Arrhenius laws 
(4). From these fits, we obtain the activation energies 
Q = 0.49 cV for D Ah Q = 0.48 cV for L> Ni , and Q = 
0.51 eV for L. The prefactors D are 1.05 x 10~ 7 m 2 /s, 
1.15 x 10" 7 m 2 /s, and 0.91 x 10" 7 m 2 /s for D Ah D m , and 
L, respectively. These values for the prefactors can be 
compared to those proposed by Eq. (4) . With the reason- 
able choices A = 3 A and v = 6THz, similar values for Do 
as in the fits are obtained, i.e. Do w 10~ 7 . 

Moreover, the expression (3) for the kinetic coefficient 
fcwF does not predict the order of magnitude correctly. To 
see this, we can compute the value of fcwF at T m « 1520 K 
using the results from the simulation. With I = 0.23 eV, 
D a w 2 x 10~ 9 m 2 /s, A 2 /6 w 1.5 A 2 , d = 3 A and / = 1, 
the value fcwF ~ 0.05m/s/K is yielded which is about one 
order of magnitude higher than the value for k, as obtained 
in our simulation. The result for fcwF is even worse if 
we replace the self-diffusion constant by the intcrdiffusion 
constant in our estimate. In this case, we find fcwF ~ 
0.25m/s/K. 

But why does the Wilson-Frenkel theory overestimate 
the speed of crystal growth? To address this question we 
propose the following scenario: We assume that the speed 
of crystal growth is limited by the atoms in the liquid- 
crystal interface region and not by the atoms in the liquid 
region where the liquid behaves like a bulk liquid. If this 
is true, one has to study the diffusion dynamics in the in- 
terface region: if diffusion in the interface region is much 
slower than in the bulk liquid, a failure of the Wilson- 
Frenkel model would be plausible since this model only 



takes into account diffusive transport of the bulk liquid. 
To check this scenario, we have simulated inhomogeneous 
systems with two crystal-liquid interfaces at the melting 
temperature T m = 1520 K. From these runs, we deter- 
mined the diffusion profiles D z ^ a (z) (a = Ni,Al) along 
the z direction that are shown in Fig. 6 together with the 
number density profile. D z , a (z) was computed from the 
long-time limit of the mean squared displacement in z di- 
rection, 



1 Ns 

D z , a (z s ) = lim — V 

t^oo iV, * — ' 



{(zi.it) - z ls (0)) 2 ) 
2t 



(5) 



where Zi s is the z coordinate of a tagged particle that was 
at time t = in one of 30 slabs that we introduced along 
the z direction, each slab having a thickness of about 2.4 A. 
N s is the number of particles in slab s (s = {1, ...,30}). 
As can be seen in Fig. 6, the interface region extends over 
5-6 atomic layers. Within this region the self-diffusion 
constants decrease roughly by about one order of mag- 
nitude. When one considers crystal growth, this slowing 
down of diffusion has to be taken into account, since the 
formation of new crystalline layers occurs in the interface 
region. This can be the reason why the Wilson-Frenkel 
model overestimates the speed of crystal growth. 

Conclusions. — Extensive MD simulations have been 
used to investigate the crystallization kinetics as well as 
the diffusion dynamics of AI50N150. Although crystal 
growth is relatively slow in this system, the simulation 
yields accurate estimates of the melting temperature and 
the kinetic growth coefficient [for the (100) orientation of 
the intermetallic B2 phase]. The small value of the lat- 
ter quantity k = 0.0025m/s/K, reveals that the growth 
kinetics of the intermetallic B2 phase is controled by dif- 
fusive mass transport. However, the classical model for 
diffusion-limited growth due to Wilson and Frcnkel docs 
not give an accurate description. We argue that this is due 
to the neglect of diffusive transport in the crystal-liquid 
interface region. Microscopic theories of crystal growth 
shall take into account the latter diffusive processes. 
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